The free surface of superfiuid 4 He at zero temperature 
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The structure and energetics of the free surface of superfiuid 4 He are studied using the diffusion 
Monte Carlo method. Extending a previous calculation by Valles and Schmidt, which used the 
Green's function Monte Carlo method, we study the surface of liquid 4 He within a slab geometry 
using a larger number of particles in the slab and an updated interatomic potential. The surface 
tension is accurately estimated from the energy of slabs of increasing surface density and its value 
is close to one of the two existing experimental values. Results for the density profiles allow for the 
calculation of the surface width which shows an overall agreement with recent experimental data. 
The dependence on the transverse direction to the surface of other properties such as the two-body 
radial distribution function, structure factor, and one-body density matrix is also studied. The 
condensate fraction, extracted from the asymptotic behavior of the one-body density matrix, shows 
an unambiguous enhancement when approaching the surface. 

PACS numbers: 67.40.Rp, 68.03.Cd, 02.70.Ss 



I. INTRODUCTION 

The interest of quantum many-body theory in inhomo- 
geneous superfiuid 4 He has been active for many years 
and continuously enriched by the achievement of new 
experimental realizations of confining geometries. 1 The 
most classical, and consequently best known, correspond 
to 4 He films adsorbed on different substrates^ and to 
4 He clusters produced by free jet expansions through very 
thin nozzles^ More recently, and still with some open 
questions, interest has been devoted to 4 He in porous 
media such as vycor or aerogel^* 6 ^ and 4 He adsorbed in 
carbon nanotubc bundlesSi^iiSiiiiiS forming quasi-one di- 
mensional structures. In common with other fluids cur- 
rently under study, the effects of a confined geometry 
upon its microscopic properties are most interesting. Ad- 
ditionally, liquid 4 He poses an unavoidable quantum be- 
havior which manifests in its own existence as a liquid 
even at zero temperature, with a superfiuid phase and 
a partially occupied Bose-Einstein state. Moreover, its 
extreme purity allows for a much cleaner extraction of 
data which is usually impossible with other liquids. 

One of the fundamental features raised by the inho- 
mogeneous situations described above is the emergence 
of the free surface of superfiuid 4 He.— For a long time, 
the main concern has been to understand how distinctive 
properties of this quantum liquid, like its superfiuid den- 
sity and condensate fraction, behave when the density 
goes to zero through a free surfaced Also the surface 
tension, the form of its density profile, the value of the 
surface width, and the role of the surface excitations in 
its dynamics have been theoretically and experimentally 
studiedAiSiiSiiLi& The most direct information on the 
density profile and interfacial width has been obtained 



by x-ray reflectivit; 
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and ellipsometrioaa measure- 



ments. The latter was carried out by Osborne^ who 
assuming a Fermi function for the profile, extracted a 
10-90 % width of 9.4 A at 1.8 K. The most recent data 



correspond to the x-ray measurement of Penanen et cUJL 
who reported a value 6.5 A at 0.45 K, a lower tempera- 
ture than the first measure performed by Lurio et alm& 
at 1.13 K who reported a wider surface of 9.1 A. The 
relative spread of the experimental results on the micro- 
scopic characteristics of the density profile points out the 
difficulties which experimental setups face to carry out 
a clean extraction of the data and its reduction to the 
ground state, i.e., to zero temperature. 

The difficulty of the experimental work on the study 
of the 4 He free surface is also reflected in the long-time 
effort to measure accurately the surface tension. Since 
the first work by Urk et al~ in 1925, the number of 
groups and different techniques used for this measure is 
unusually large. Recent data starts with lino et al. 2i 
who measured the surface tension using the surface-wave 
resonance method obtaining an estimation at zero tem- 
perature of a = 0.257 HA -2 . Later on, Roche et alm^ 
obtained a slightly larger value a = 0.274 KA~ 2 from 
the frequencies of capillary waves with wavelength in the 
micron range. In this second work, the difference in the 
values obtained for a is imputed to a possible inaccuracy 
in the work of lino et alM> related to the treatment of 
the meniscus at the edge of the box. However, this argu- 
ment was refuted posteriorly by Nakanishi and Suzuki 26 
who confirmed the result obtained in Ref. Recently, 
Vicente et alSL performed a new measure of the surface 
tension using the vibration modes of levitated 4 He drops, 
and the result was in perfect agreement with the larger 
value obtained by Roche et al.S^ 

The surface tension, surface width, and density pro- 
file of the free 4 He surface have been calculated us- 
ing both density functionaliSiS 2 ^ and microscopic the- 
or y > i5 i 30 i 3i 4 32 J 33 Density functional theory relies on func- 
tional forms of the energy which depend on the one-body 
density and some parameters which are adjusted to repro- 
duce selected experimental data. These functionals have 
been progressively improved by the inclusion of non-local 
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contributions and specific terms to reproduce the experi- 
mental static response function™ Contrarily to the first 
models, in which the surface tension was considered an 
input to fix parameters, the most modern ones are able 
to give predictions for a. The results for the surface 
tension are in an overall agreement with the experimen- 
tal data from Ref. |2^, and the surface width is ~ 6 A. 
The methodology of microscopic approaches is different 
since the starting point in all of them is the Hamilto- 
nian of the liquid containing realistic interparticle inter- 
actions. In this respect, the variational method has been 
specially useful due to the hard-core of the He-He inter- 
action. In fact, the most extensive work in the field of in- 
homogeneous 4 He has been carried out in the variational 
framework of the hypernetted chain (HNC) and corre- 
lated basis function (CBF) theories. Of special relevance 
in this field is the continued work by Krotscheck and 
co-workersj24 who have carried out a nearly exhaustive 
study of 4 He in reduced geometries both for the ground 
and the excited states. 

Quantum Monte Carlo method™^ have also been ap- 
plied to the study of the free surface of liquid 4 He. Us- 
ing a slab geometry, Valles and Schmidt™ calculated 
the ground-state properties of the surface by means of 
the variational Monte Carlo (VMC) method with a trial 
wave function containing two- and three-body correla- 
tions. They extracted a surface tension a = 0.272 KA -2 
(again very close to data from Ref. 25) from an ana- 
lytical fit to the dependence of the energy with the in- 
verse of the surface density. Galli and Reatto^i studied 
a 4 He slab using a shadow wavefunction and the VMC 
method; they obtained a = 0.31(1) KA -2 and a sur- 
face thickness ~ 5 A. The variational constraints can 
be removed by applying the essentially exact diffusion 
Monte Carlo (DMC) and Green's function Monte Carlo 
(GFMC) methods. At present, the only application of 
these methods to a slab corresponds to a GFMC calcu- 
lation by Valles and Schmidt who estimated a surface 
tension of 0.265 KA" 2 and a small value for the width 
(~ 4 A) probably due to the small number of particles 
used in the simulation. Recently, Draeger and CeperleyS 
carried out a similar study but at finite temperature us- 
ing the path integral Monte Carlo (PIMC) method. This 
work does not contain results for a and surface width 
since its main motivation was the study of the enhance- 
ment of the Bose-Einstein condensate fraction at the sur- 
face. 

In the last years, a renewed interest in the study of 
the 4 He surface has emerged since the first theoretical 
observation by Lewart et aZ™ of a large increase of the 
condensed fraction near the surface of He drops. The 
relevance of this result in the long way for searching a 
Bose-Einstein condensed state was put forward by Griffin 
and Stringar™ who studied the low density surface using 
a generalized Gross-Pitaevskii equation. A subsequent 
work by Galli and Reatto^ using the VMC method with 
a trial wave function based on the shadow model, pointed 
out that the introduction of fluctuations on the surface, 



mainly the zero point motion of ripplons, can reduce sig- 
nificantly the enhancement of no. Their results show an 
increase of no up to a maximum value 0.5 and a de- 
crease to zero when the density p(r) approaches zero. 
In order to shed light on these two different predictions, 
Draeger and Ceperley 3 ^ calculated no in a 4 He slab us- 
ing PIMC. Their results show also a clear enhancement 
of the condensate fraction up to values larger than 0.9 
and a small decrease when approaching the outer part of 
the surface. Therefore, fluctuations can arise and reduce 
no but their effects could be si gnifi cantly smaller than 
in the VMC calculation of Ref. |3(J A conclusive state- 
ment would require from a direct experimental measure 
which, for example, in the bulk has not yet been possible. 
Nevertheless, at present there are two experiments that 
seem to confirm the increase of no. The first one was 
derived from a quantum evaporation experiment^ and 
the second and more recent, form deep-inelastic neutron- 
scattering on thick layers of 4 He on an MgO substrate™ 

In the present paper, we present a DMC calculation of 
the ground-state properties of a free 4 He surface at zero 
temperature. To some extent, our work represents an up- 
date of the GFMC calculation of Valles and Schmidt 32 
in the sense of considering a more accurate interatomic 
potential and a larger number of particles in the simula- 
tion to achieve thicker slabs and therefore a more realistic 
study of the surface. Moreover, we have included in the 
present study topics like the density profile form, two- 
body distribution functions and structure factors, and 
mainly the study of the one-body density matrix which 
were not analyzed in Ref. lU The exact character of 
the DMC method, within the statistical noise, allows for 
an accurate study of all the surface properties without 
the bounds imposed by a variational treatment. The 
methodology used in the analysis is, regarding the main 
inputs, the same that we have used previously in the 
study of bulk liquid 4 He and 3 He, and that has shown an 
overall agreement between theory and experiment™ 

Our paper is organized as follows. In the next section, 
we review the DMC method which works with a second- 
order approximation for the short-time Green function. 
The specific details of the simulation for the slab geom- 
etry are also discussed, with special emphasis on the re- 
quirements for the trial wave function used for impor- 
tance sampling. The results obtained, for both structure 
and energetics of the slab, are presented in Section III. 
Finally, Section IV comprises additional discussion on the 
results and the main conclusions. 



II. DIFFUSION MONTE CARLO AND THE 
SLAB GEOMETRY 

The DMC method is nowadays a well established tool 
for solving the many-body imaginary-time Schrodinger 
equation, 
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where R = (rx,...,rjv) is a 3A-dimensional vector 
(walker) and t is the imaginary time measured in units 
of h. The time-dependent wave function of the system 
^(R, t) can be expanded in terms of a complete set of 
eigenfunctions 4>i (R) of the Hamiltonian, 



*(R,t) = ^c n e X p[-{Ei-E)t] &(R) 



(2) 



where Ei is the eigenvalue associated to <f>i(R). The 
asymptotic solution of Eq. JQ), for any value E close 
to the energy of the ground state and for long times 
(t — > oo), gives (f>o (R), provided that there is a nonzero 
overlap between ^(R, t — 0) and the ground-state wave 
function (f>o(R). 

An efficient solution of Eq. requires from the use 
of importance sampling. It is introduced by solving the 
Schrodinger equation for the wave function 



/(R,t) = ^(R)*(R,t) , 



(3) 



ip(R) being a time-independent trial wave function able 
to describe the system at a variational level. Introducing 
the Hamiltonian of the system 



H = 



2m 



V 



R 



V(R) 



(4) 



Eq. 0} turns out to be 
df(R,t) 



Of 



= -DV 2 n f(R,t)+DV R (F(R)f(R,t)[ 



(E L (R)-E) f(R,t) , 



(5) 



with D = h 2 /(2m), E L (R) = ^(R^H^R) the local 
energy, and 



F(R) = 2^(R)- 1 V R V'(R) 



(6) 



is the drift force which guides the diffusion process. 

The r.h.s. of Eq. (JSJ) may be written as the action of 
three operators Aj acting on the wave function /(R, t), 



g/£M) 

dt 



= (A 1 + A 2 + A 3 )f(R,t) = Af(R,t) (7) 



The first one, Ai, corresponds to a free diffusion with 
a a diffusion coefficient D; A 2 acts as a driving force 
due to an external potential, and finally A 3 looks like 
a birth/death term. Equation 0) is transformed to the 
integral form 

f(R', t + At) = J G(R', R, At) f(R, t) dR , (8) 

by introducing the Green function 

G(R, R, At) = (R | cxp(-AAt) \ R) . (9) 

with A = Ax + A 2 + A 3 . 

In the DMC method, Eq. (0 is iterated repeatedly 
until to reach the asymptotic regime f(R, t — > oo), a limit 



in which one is effectively sampling the ground state. In 
our implementation of the method we use a second-order 
expansion (quadratic DMC) ^ 

exp {-AAt) = (10) 

exp^-A 3 ^ exp^-A 2 ^ exp(-AiAt) 

x exp (-A 2 ^£) exp (-A 3 ^f) + O ((At) 3 ) , 



which has proved to be a good compromise between al- 
gorithmic complexity and efficiency. 

The study of the free surface is made by simulating 
a slab which grows symmetrically in the z direction and 
with periodic boundary conditions in the x— y plane. The 
slab geometry is probably the best suited one to analyze 
the free surface since it does not present neither the sub- 
strate influence of films nor the curvature effects of drops. 
It is also true that the existence of two surfaces can in- 
troduce some residual influence between them. However, 
we have checked in our simulations that this effect, at 
least for the larger slabs, is completely negligible. 

The implementation of DMC requires of a model for 
the trial wave function ip(R) with the basic physical fea- 
tures of the system under study. In the present case, one 
a priori knows that ^(R) becomes zero in two situations: 
when two atoms get closer than the core of its interaction, 
and when an atom tries to escape from the surface. No- 
tice that at zero temperature no vapor is present out of 
the liquid. Having in mind these arguments, we consider 



^(R)=Vj(R) 0(R) 



(11) 



with a Jastrow correlation factor accounting for dynam- 
ical correlations 



N 



^(R)=n/(^o » 



(12) 



and a confining term, factorized in the form 

N 



(13) 



i=l 



The two-body correlation function is the same we used 
in the study of the bulk liquid^ It was proposed by 
Reattoa^ and incorporates, in an approximate way, the 
medium range behavior of f(r) observed in the Euler- 
Lagrange functional optimization, 



f(r) = exp 



- (- 

2 \r 



L 

2 exp 



r — A 
A 



The confining function is chosen of Fermi type, like in 
Refs. andll 



h(z) = {1 + exp[fc( \z - z cm | - z )]} 1 , 



(15) 
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with parameters k and zq controlling the width and loca- 
tion oi the interlace, respectively. Possible and spurious 
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particle coordinate z in Eq. (JTSJ. 
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1.1161 


22.10 


-6.766(11) 


13.075(40) 


III. RESULTS 




oo 




-7.267(13) 


14.32(5) 



The properties of the free surface of liquid 4 He have 
been studied by performing DMC calculations with an 
increasing number of particles JV=54, 108, 162, 216, 324, 
and a fixed x — y area of A — 290.3 A 2 . The interatomic 
potential is the HFD-B(HE) model proposed by Aziz et 
alm^ that has proven its high accuracy in previous DMC 
calculations of the density dependence of the pressure 
and speed of sound in bulk 4 He and 3 He liquids^ As 
a matter of comparison with the GFMC calculation of 
Valles and Schmidt, 32 and in order to establish the influ- 
ence of the potential in the surface properties, we have 
calculated the energy of the system and the surface ten- 
sion using the HFDHE2 model, also proposed by Aziz 
and collaborators^ 

In order to take into account size effects, due to the fi- 
nite value of the simulation area A and periodic boundary 
conditions in the x—y plane, we have added to the energy 
tail corrections which are externally calculated assuming 
that the particles are uncorrelated for distances in this 
plane larger than half of the simulation square. The an- 
alytical expressions for the tails are reported in Ref. l32l 
Applying this analysis, we have verified that the energy 
per particle when doubling simultaneously A and N does 
not change within the current statistical noise. Possible 
bias in the DMC calculations coming from the time step 
and the mean population of walkers are also under con- 
trol working, in both cases, with values which are well in 
the asymptotic regime. 

A final key point of a DMC simulation concerns the val- 
ues of the parameters in the trial wave function (|14I15J) . 
This has been dealt by performing a series of VMC calcu- 
lations and searching the optimal set. As expected, the 
only parameter which clearly shows a dependence on the 
number of particles, i.e., the size of the slab, is zq\ the 
particular values are reported in Table I. The rest of the 
set has shown negligible dependence on N; their values 
are b = 3.067 A, L = 0.2, A = 5.112 A, A = 1.534 A, and 
k = 1 A" 1 . 



A. Energies and surface tension 

The DMC results of the energy per particle as a func- 
tion of the number of particles included in the simulation, 
with a fixed area A, is reported in Table I. The coverage 
densities n c , defined as 



TABLE I: Total (E/N) and kinetic (T/N) energies per parti- 
cle as a function of the number of particles N and the cover- 
age density n c . E 1 /N is calculated with the old Aziz poten- 
tial (HFDHE2)(— the results in squared parenthesis are taken 
from Ref. I32I Figures in parenthesis are the statistical errors. 
20 is a variational parameter present in the correlation factor 
h(z) lO. The bulk result is taken from Ref. 



with p(z) the density profile, are also contained in the Ta- 
ble. When n c increases, the energies approach the bulk 
valued which appears in the last row to help in the com- 
parison. In the last column, we have included results for 
the energy using the old Aziz (HFDHE2) potential. As it 
is known from previous DMC calculations for the bulk, 42 
the energies obtained with the HFDHE2 potential are 
slightly smaller (in absolute value), the difference with 
the HFD-B(HE) model being mainly due to the increase 
in the potential energy. Results for the kinetic energy, 
which are the same for the two potentials within error 
bars, are included in the Table showing convergence to 
the bulk value (14.32 K)i£ Enclosed in squared brackets 
in the Table there are GFMC results from Ref. |3^ cal- 
culated at the same coverages; the agreement with our 
results is better for the larger coverage case but the lack 
of data at larger n c leaves the comparison rather incom- 
plete in the most interesting region, i.e., for the largest 
slabs. 

The energies per particle, as a function of the inverse 
of the surface coverage (l/n c ), are shown in Fig. 1. As 
a limiting case (n c — oo), the energies corresponding to 
the bulk for the two potentials are also plotted. Under 
general assumptions it is expected that the energy per 
particle in the slab follows a polynomial expansion in 
terms of l/n Cl 29 



2a 

lie 



O 



1 



(17) 



N 
A 



dz p(z) 



(16) 



In this expansion, a is the surface tension and the pa- 
rameter 7 has been related to the long-range interaction 
term of an hypothetical 4 He substrate^ In Fig. 1 (top), 
the energies for slabs with different number of atoms are 
shown for values l/n c < 3 A 2 . In this regime, corre- 
sponding to the largest slabs, the energy is a linear func- 
tion of l/n c . A linear fit to the data, for both inter- 
atomic potentials, gives x 2 /v = 1 — 1.5. From the linear 
coefficient we estimate the surface tension; for the most 
accurate HFD-B(HE) potential a = 0.281(3) KA~ 2 , and 
a = 0.278(3) KA~ 2 for the old Aziz potential. Within 
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FIG. 1: Energies per particle as a function of the inverse of 
the surface coverage. Circles and squares stand for the cal- 
culations with the HFD-B(HE) and the HFDHE2 potentials, 
respectively. Top: Large n c regime; solid and dashed lines 
correspond to numerical fits HI 711 with 7 = 0. Bottom: All 
the coverages; the lines are numerical fits l|17|l to the DMC 
data with 7 7^ 0. 



the statistical errors, the surface tension is the same for 
both potentials implying that changing the potential in- 
troduces a shift in the energy which is to a large extent 
independent of the number of atoms. In the fit, the DMC 
energy of the bulk phase, corresponding to l/n c = 0, has 
been included. However, it turns out that the inclusion 
or not of the bulk energy in the fit does not modify the 
result for the surface tension, and that the extrapolated 
energy to l/n c = matches the bulk result, a feature 
that gives us additional confidence on the accuracy of 
the slab energies. 

In order to describe the whole data set it is necessary 
to introduce in the fit higher-order terms (|1 7|) . To this 
end, we have introduced a next-order term with 1 /n 2 or 
1/njj, and tried to discriminate which of the two options 
is more favored by our data. Using the third-power, as 
suggested by Szybisz 29 ifTTjl . the resulting value for x 2 l v 
is 2.5, a value significantly smaller than the one achieved 
with a quadratic law (8.0). Therefore, and in spite of not 
being completely conclusive, our data support Szybisz's 
analysis. The corresponding fits are shown together with 
the DMC data in Fig. 1 (bottom). The surface tension 
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FIG. 2: Density profiles for an increasing number of particles 
in the slab. From left to right, they correspond to N =80, 
108, 162, 216, and 324 atoms. The dashed line shows the 
equilibrium density of the homogeneous liquid, po = 0.02186 

A" 3 . 



from this fit is a = 0.291(4) KA 2 , nearly compatible 
with our best prediction from the linear fit to the largest 
n c values. On the other hand, 7 = -0.0023(2) KA" 6 
roughly half the approximated value obtained by Szybisz. 

The present result for the surface tension is a bit larger 
than the one obtained by Valles and Schmidt in a GFMC 
calculation with the old Aziz potential (HFDHE2), a — 
0.265(6) KA _2 iS That difference can be attributed to 
the fact that we have available data corresponding to 
large n c values, which is crucial for the estimation of er, 
and somehow also to the use they made of a second- 
order polynomial fit. Chin and Krotscheck2& obtained 
a = 0.284 kA~ 2 in an Eulcr-Lagrange optimized HNC 
calculation of 4 He clusters. Density functional theory in- 
cluding non-local terms are predictive for the surface ten- 
sion^ the different functionals in the literature 2 ^ show 
results for the surface tension in the range 0.272-0.287 

kA- 2 . 



B. Density profile 

Results for the density profiles of the slab, with differ- 
ent number of atoms in the simulation cell, are shown in 
Fig. 2. In order to eliminate any residual bias, due to 
the importance sampling trial wave function, pure esti- 
mators have been used^I When N increases, the surface 
moves along z and the density in the central part in- 
creases progressively towards the equilibrium density of 
bulk 4 He (p = 0.02186 A" 3 ). This density is nearly 
reached for our largest slab corresponding to JV = 324 
atoms. The oscillation in density which appears in the 
inner part of the slab does not have any physical mean- 
ing since its amplitude is compatible with the statisti- 
cal noise there. The existence of stable oscillations in 
the surface has been discussed in several works after the 
first proposal by Reggei 4 ^ Within the present accuracy, 
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FIG. 3: Symetric vs. asymetric fit to the 4 He density profile. 
Points with errorbars are the computed values for the slab 
with N = 324. The dashed and solid lines correspond to a 
symetric (8 = 1 in Eq. JTHJ) and an asymetric model (5 ^ 1), 
respectively. 



our results exclude this possibility. Nevertheless, the re- 
sults contained in Fig. 2 show a slight shoulder in the 
inner part of the surface (hardly observable for the small- 
est slabs) that could resemble a Regge oscillation. This 
property of the He density profile has been observed 
previously in semi-infinite matter and clusters^&Sl 

The different slope of the density profile in the inner 
and outer parts of the surface has been observed and 
discussed in previous works. This feature is directly re- 
lated to the possible asymmetry of JLlSiSft In or- 
der to guide our analysis on this point, and trying to be 
quantitatively accurate, we have fitted to the DMC data 
functions of the type 



po 



(1 + e /S(«-*o)) < 



(18) 



If S = 1 in Eq. I |18|l the density profile is symmetric; oth- 
erwise, the model is asymmetric with a degree of asym- 
metry governed by S. Using \ 2 as a goodness parameter, 
we have observed that the DMC data for the larger slabs 
(N > 162) is best reproduced by an asymmetric fit. On 
the contrary for the smaller and less realistic slabs, the 
quality of the fit is not improved if 8 ^ 1. In Fig. 3, sym- 
metric and asymmetric fits to the DMC data on the slab 
with TV = 324 are shown on top of the simulation results. 
The differences between the two fits are certainly small 
but the asymmetric model matches better the DMC re- 
sults, especially in the inner part of the surface where 
the slope of the profile is less pronounced. An additional 
argument of consistency towards the assessment of an 
asymmetric profile is the stability in the value of S which 
is the same for the three larger slabs, S = 1.91(15). 

One of the most relevant properties characterizing the 
surface is its width, usually measured as the length (W) 
over which the density decreases from 90 to 10 % of the 
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0177(21 
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7.75(5) 
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4.01(5) 


10.09(5) 


162 


0.0198(2) 


4.57(5) 


14.19(5) 


216 


0.0206(2) 


4.76(5) 


18.06(5) 
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0.0217(2) 


5.07(5) 


25.74(5) 



TABLE II: Central density (po), surface width (W), and size 
of the slab measured with Z1/2, as a function of the number 
of atoms N. 




FIG. 4: Interfacial width W as a function of the size of the 
slab, measured in terms of Z\ n (solid circles). The open circles 
correspond to experimental data from Penanen et al.w~ The 
solid line is the fit 1191 considering only our theoretical data 
and the dashed line corresponds to the same analytical form 
but including also the experimental points. 



inner constant value. Theoretical calculations of the 4 He 
free surface at zero temperature using microscopic the- 
ory 15 i 30 i 31 i 32 i 33 or density functional approacheo 16 i 28 i 29 
predict values in the range 5-7 A. These results are in 
overall agreement with the most recent experimental data 
at T — 0.45 K, which point to a width of 6.5 A for thick 
(125 A) filmsi^i Our results for the width are reported in 
Table II as a function of the number of particles in the 
slab N; the size of each slab is given in terms of 21/2, 
defined as the value of z in which the density is half the 
density in the inner part po . The results contained in the 
Table have been extracted from the DMC density pro- 
file by performing a fit with the asymmetric form given 
by Eq. (JTSJ. If the fit is symmetric (5 = 1), W is re- 
duced with respect to the values reported in the Table 
but the differences are only ~ 0.1 A since the degree of 
asymmetry is rather small. 

The main prediction of the present study would be 
the value of W corresponding to a free surface, with 
a central density equal to the bulk equilibrium value 
Po ulk = 0.02186 A~ 3 . However, the values reported in 
Table II show a dependence on the size of the slab, even 
for the largest ones, which makes difficult that estima- 
tion. In order to shed light on this dependence, and try- 
ing to analyze possible extrapolations to larger slabs, the 
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FIG. 5: 2-dependence of the total (E/N), kinetic (T/N), and 
potential (V/N) energies in the slab geometry. The behavior 
of the energies is compared with the density profile (right 
scale). The data corresponds to a calculation with N = 108 
atoms. 




FIG. 6: Two-body radial distribution function for different 
values of z, g(ru, z). Going from bottom to top in the height 
of the main peak the data corresponds to z — 11, 9.5, and 2 
A for a slab with 324 atoms. 



results of the width W are shown as a function of zi/ 2 in 
Fig. 4. After testing different analytical forms, we have 
arrived to a rather simple model that describes quite ac- 
curately the DMC data, 



is proportional to the analytic expressionist 

E/N(z) = T/N(z) 



lim - 1 V2 ^P ■ (20) 



2m ^/p{z) 



W(z) =W 00 (l- e" 6 ^) 



(19) 



Woo being the predicted width of the free surface. The 
value obtained is Woo = 6.3(4) A. In Fig. 4, experimental 
results from Penanen et ai.j2i corresponding to thick 4 He 
films adsorbed on a solid substrate at T — 0.45 K, are 
also shown. If the DMC and experimental data are fitted 
altogether, according to the empirical law JTJJ, the inter- 
facial width of the free surface becomes Woo = 6.2(1) A 
which is compatible with the extrapolated value consid- 
ering only the theoretical results. Both extrapolations 
are in agreement with the width measured directly for 
the thickest film of 125 A, W = 6.5 ± 0.5 A. 

The transition from a constant density in the central 
part of the slab to a zero density in the outer part of the 
surface implies a z-dependence in other magnitudes, in 
particular in the total and partial energies. In Fig. 5, 
this z dependence is shown in comparison with p(z) for 
a N = 108 slab. The three functions E/N(z), V/N(z), 
and T/N(z) are calculated by averaging the energies of 
particles located between z and z + Az according to a 
uniform grid with Az ~ 0.1 A and, therefore, they are 
not exact estimations. Nevertheless, the residual bias 
coming from the trial wave function is expected to be 
small and the qualitative behavior would not change from 
the one reported in Fig. 5. As one can see, the potential 
energy evolves from a constant value in the inner part 
(corresponding to the bulk value at the central density 
of this slab) to zero when p(z) — » 0. Also the total and 
kinetic energies in the central part of the slab correspond 
to the bulk values but, approaching the surface, the total 
and kinetic energies reach the same constant value which 



C. Distribution and structure functions 

The slab geometry provides the opportunity of study- 
ing the density dependence of the spatial structure in 
a quantum liquid from bulk to zero densities. The z- 
dependence of the structure functions is analyzed by per- 
forming slices of variable width Az (larger in the center 
of the slab and smaller in the surface) in which the den- 
sity can be considered like a constant. The structure of 
the fluid is then calculated as a function of z, accumulat- 
ing data corresponding to atoms localized between z and 
z + Az. Following this scheme, we have calculated the 
two-body radial distribution function g(r\\, z) and the re- 
sults obtained, for a selected set of z values, are reported 
in Fig. 6. The functions g(ru,z) follow the evolution of 
the density: by increasing z (decreasing p) the height of 
the main peak decreases and moves to larger rii values, 
reaching in the outer part of the surface a very low den- 
sity regime where g(r\\, z) resembles the typical result for 
a dilute gas. 

Additional information on the structure of the sys- 
tem, and on its excitation modes, are contained in the 
z-dependent static structure factor. In its general form 
it is defined as 



S(k\\,z,z') = 



[N(z)N{z')} 1/2 



(Ph n (z) P-k n (z')) , 



(21) 

Phi] ( z ) — Si=i e lfc ii' ri being the fluctuation density op- 
erator for all the particles N(z) in the slice (z, z + Az), 
and fen a momentum in the x — y plane. Our calculations 
have addressed only the diagonal part of this function, 
z = z', which in a simplified notation are called S(k\\, z) 
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FIG. 7: Static structure function for different values of z, 
S(k\\, z). Going from bottom to top in the height of the main 

peak the data corresponds to z = 11, 9.5, and 2 A for a slab FIG. 8: One-body density matrix in a 4 He slab with N = 324. 
with 324 atoms. 



in this work. Results for S(k\\,z), corresponding to the 
same slices chosen for g(ru,z) in Fig. 6, are shown in 
Fig. 7. The evolution of S(k\\, z) with z follows the same 
trends observed in Fig. 6 for g(ru,z), i.e., a behavior 
essentially determined by the local density in each slice 
p(z). The behavior when k\\ — > is known to be dif- 
ferent in the central part and in the surface of the slab. 
In the central part, where the mean density is constant, 
the static structure factor must show a linear behavior 
with fc| | , corresponding to the phonon mode of the homo- 
geneous liquid. On the contrary, in the surface S(k\\,z) 
diverges as \/^/k\\ due to the presence of ripplons. 51 The 
DMC results shown in Fig. 7 do not show signatures of 
ripplons due probably to restrictions imposed on the low- 
est fc|| value accessible with a finite number of particles. 
Recent results for the structure function^ using path 
integral Monte Carlo (PIMC) with a number of particles 
larger than the one used in this work, show a change of 
behavior for the smallest fc|| values and for the most ex- 
ternal slices that could be compatible with the singular 
behavior induced by ripplons. Notwithstanding, the size 
of this signal observed in the PIMC results is appreciably 
reduced with respect to a VMC calculation, using shadow 
wave functions, which emphasized the importance of rip- 
plons for a correct description of the 4 He surface. 36 



D. Condensate fraction 

The presence of off-diagonal long-range order 
(ODLRO) in a Bose liquid is one of its more essential 
properties^ The measure of ODLRO in the system 
is quantified by the Bose-Einstein condensate fraction 
no, i.e., the fraction of particles occupying the zero- 
momentum state. This information is contained in the 
one-body density matrix 



N J dr 2 ...r N *(ri,r 2 , . . . ,r N )^(r[,r 2 , . ■ .,r N ) . 

The geometry of the slab makes possible to express the 
one-body density matrix as a function of rn, z, and z', 
p{r\\, z, z'). Restricting the analysis to the particular case 
z = z' , the long range behavior of p(r\\, z) gives the local 
condensate fraction no(z), 



n (z) 



lim 



P(r\\,z) 
p(Q,z) 



(23) 



The function p(r\\, z) is sampled in DMC following the 
standard procedure, i.e. through the expected value 



ip(r[,r 2 , ...,r N ) 
tp(ri,r 2 ,...,r N ) 



(24) 



p(ri,r[) 



(22) 



accumulating data for a set of slices in the z direction as 
explained in the previous subsection. Results for p{r\\, z), 
corresponding to the largest slab (N — 324), are shown 
in Fig. 8. In this 3D picture, one can see the constant 
regime achieved for any z value when ru — » oo, which 
corresponds to the condensate fraction no(z) l|23|l . Fol- 
lowing the evolution of tiq(z) with increasing z, one ob- 
serves a large enhancement approaching the surface up 
to values close to 1 which would correspond to a fully 
Bose-condensed state. 

The change in the condensate fraction as the surface is 
approached is more clearly shown in Fig. 9. In the inner 
part of the slab, where the mean density is very close 
to the bulk equilibrium density, no(z) is constant and 
equal to 0.095(5). When the density starts to decrease, 
no(z) increases up to almost 1. This result is not sur- 
prising since it is well known from bulk calculations that 
no decreases monotonically with the density. Possible 
density fluctuations in the surface, due to the presence 
of ripplons (surface excited modes), that can suppress 
this nearly Bose-condensate state are not present in our 
calculations because the goal of the present work is the 
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studied could be rather well described by the law 
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0.000 



z(A) 



FIG. 9: z dependence of the condensate fraction in a 4 He slab 
with N — 324 (solid line, left scale). The condensate density 
po(z) (long-dashed line) and density profile p(z) (short-dashed 
line) are also shown (right scale). 




0.4 0.6 

P(z)/P 



FIG. 10: Condensate amplitude \fna\p{z)\ as a function of 
the density p(z)/po, with po the bulk equilibrium density. The 
points correspond to estimations of the condensate fraction 
for all the slabs studied. The solid line corresponds to the 
analytical fit 125H to data in the density range 0.1 < p(z)/po < 
1. 



ground state of the free 4 He surface. At low tempera- 
tures, the influence of ripplons in tiq(z) has been studied 
by Draeger and Ceperley using the PIMC method^ their 
results show an at most 20 % reduction of no, with ef- 
fects only at extremely low density where, on the other 
hand, the statistical signal is rather poor. Figure 9 also 
contains the evolution with z of the condensate density, 
Po{ z ) — n o{ z )pi. z )i which shows a clear maximum at 
approximately z\/2 and then becomes equal to the to- 
tal density when p(z) — > 0. The resulting picture is in 
an overall agreement with the VMC results of Lewart 
and Pandharipande^i and with recent DMC studies of 
trapped hard-core bosons by DuBois and GlydeS 

Collecting data for hq as a function of the density, and 
for the different slabs calculated, one can try to elucidate 
if they can be well reproduced by a single analytical form. 
This was already analyzed by Lewart et alM> in a VMC 
calculation of natural orbitals in 4 He clusters. They ob- 
served that the condensate fraction for all the clusters 



y/n {p(z)} = A-B- 



Po 



(25) 



The function yj n§\p(z)\ can be seen as a condensate 
amplitude assuming that po(z) ~ |^E'o( 2: )| 2 ! with &o(z) 
the condensate wave function. The present results for 
y/ no[p(z)], including the data for all the slabs, are re- 
ported in Fig. 10 as a function of p(z)/p a , p a be- 
ing the equilibrium density of bulk 4 He. In the range 
0.1 < p(z)/p < 1 all the data shows the same linear be- 
havior that can be very well fitted by the empirical law 
J23> with values A = 0.935 and B = 0.636. As can be 
seen in Fig. 10, the DMC data departs from this linear 
behavior at very low densities p(z)/po < 0.1, approaching 
faster to the zero value at zero density. 



IV. DISCUSSION 

Using a slab geometry we have presented in this work 
a DMC study of the ground-state properties of the free 
4 He surface. The interatomic potential is highly accu- 
rate, and in the past has allowed for an excellent repro- 
duction of the experimental equation of state in the bulk 
phased We believe that the combination of an essen- 
tially exact method and this very well know interaction 
between helium atoms guarantees the accuracy of a mi- 
croscopic description of the kind intended here. At this 
point it is interesting to compare the main results ob- 
tained with available experimental data. 

The value of the surface tension of superfluid 4 He at 
zero temperature is one of the main results. Our pre- 
diction, a = 0.281(3) KA~ 2 , must be compared with 
the two available and different experimental measures, 
a = 0.257±0.001 HA- 2 and a = 0.274±0.002 kA" 2 ^ 4 * 2 ^ 
Both experimental determinations correspond to zero 
temperature extrapolations using measured values for 
a(T) and the Atkins ripplon law for a(T) at very low 
temperatures^ Our theoretical result is larger than both 
measures but it is nearly statistically compatible with the 
value er = 0.274 KA~ 2 obtained by Roche et al.M It is 
worth noticing that most of theoretical estimations of cr 
agree with this larger experimental value. 

A second result which can be compared with experi- 
ment is the surface width W . Obviously, the most rele- 
vant result would be the real width of the surface, i.e., 
the asymptotic width for thick slabs. In our results, 
W shows a dependence with the size of the slab and, 
even for the largest one, this asymptotic regime is not 
achieved. For the N = 324 slab the width is 5.07(5) 
A. Using a reasonable fit to the data we can estimate 
the asymptotic surface width, the value obtained being 
Woo = 6.3(4) A which is in good agreement with a recent 
experimental measure at T = 0.44 K for a 125 A film, 
W = 6.5 ± 0.5 A. 21 On the other hand, the density pro- 
files, specially the ones with a greater number of particles, 
show a slight asymmetry with a deeper slope in the outer 



10 



part of the surface. This possible asymmetry, observed 
also in other theoretical approaches jl&Si was analyzed 
in an x-ray specular reflectivity experiment by Lurio et 
alm& but uncertainties in the scattering amplitude made 
impossible a clean answer to this point. 

Last but not least there is the question of the enhance- 
ment of the condensate fraction in the surface. The ex- 
perimental measure of no in liquid 4 He has been elusive 
for long time due to its small valued Nowadays, deep 
inelastic neutron scattering has largely improved both 
its efficiency and data analysis and therefore there is 
much more confidence on the results obtained. Recently, 
Pearce et alm& have carried out neutron scattering ex- 
periments on liquid 4 Hc adsorbed in thick layers on an 
MgO substrate. The analysis of the dynamic structure 
function obtained for different number of layers shows un- 
ambiguously an increase of the Bose condensate when the 



number of layers is reduced. The data obtained are some- 
how not completely accurate at the quantitative level but 
the main conclusion of the experiment, i.e., the enhance- 
ment of no in the surface, is in agreement with our (and 
others 33,37 ) theoretical calculations. It is worth mention- 
ing that signatures of a large condensate fraction have 
been also observed by Wyatls^ in quantum evaporation 
experiments. 
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